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Abstract:  We  examine  the  use  of  a  single-fluid  model  with  the  Peng  Robinson  equation  of 
state  to  model  supercritical  and  transcritical  flows  with  combustion.  Both  non-reacting  and 
reacting  flows  are  considered  to  understand  the  modeling  challenges.  Several  modifications 
of  the  equations  of  state  and  mixture  rules  are  tested  and  shown  to  have  varying  strengths 
and  weaknesses.  No  method  works  uniformly  well  for  both  non  reacting  and  reacting  flows, 
although  the  use  of  the  Peng-Robinson  model  with  Amagat’s  mixture  rule  and  modified 
compressibility  factors  appears  to  be  robust  and  reasonably  accurate  for  supercritical  and 
transcritical  combustion. 


Nomenclature 


am  Cubic  parameter 

A  Cubic  parameter 

bm  Cubic  parameter 

B  Cubic  parameter 

kU  Interaction  parameter 
pc  Critical  pressure 

pr  Reduced  pressure  (p/pc) 

Tc  Critical  temperature 
Tr  Reduced  temperature  (T/Tc) 
Vm  Molar  volume 

Xi  Mole  fraction  of  species  i 

Z  Compressibility 

uj  Acentric  factor 


I.  Introduction 

Liquid  rocket  and  gas  turbine  engines  operate  at  high  pressures.  For  gas  turbines,  the  combustor  pressure 
can  be  60  —  100  atm,  while  for  rockets  it  exceeds  100  atm  and  can  be  as  high  as  300  —  500  atm.  Under 
these  conditions,  the  propellants  are  typically  injected  at  supercritical  pressures,  but  are  often  at  subcrit- 
ical  temperatures.  The  situation  is  further  complicated  since  the  propellants  are  injected  into  a  mixture 
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environment  comprised  of  reactants  and  products.  Whether  or  not  the  mixture  is  supercritical  (i.e.,  above 
the  vapor  dome)  is  a  function  of  the  composition  in  addition  to  the  local  pressure  and  temperature.  As  a 
supercritical  fluid  is  injected  and  mixes  with  the  ambient  fluid  it  may  become  locally  subcritical  even  if  the 
temperature  and  pressure  remain  constant.  The  subcritical  and  supercritical  regions  of  the  flow  can  have 
vastly  different  thermodynamic  and  transport  properties,  which  can  influence  mixing  and  combustion  dy¬ 
namics  in  the  combustor.  Eventually,  regardless  of  whether  the  reactants  are  in  a  subcritical  or  supercritical 
state,  combustion  occurs  and  the  resulting  high  temperatures  drive  the  system  to  a  supercritical  state.  The 
existence  of  such  a  combination  of  states  may  be  referred  to  as  a  transcritical  state.  Computational  modeling 
of  such  phenomena  is  extremely  challenging  and  is  the  subject  of  the  current  article. 

A  standard  approach  for  supercritical  flows  is  to  treat  the  multicomponent  mixture  of  species  as  a  dense 
fluid  using  a  real  gas  equation  of  state  such  as  the  Peng- Robinson  cubic  equation  of  state.1  While  the  model 
is  clearly  valid  for  flows  that  are  well  above  the  vapor  dome,  their  validity  in  the  transcritical  region,  where 
one  or  more  of  the  fluid  components  may  drive  the  system  to  a  saturated  state,  is  not  well  established. 
The  justification  for  the  use  of  the  dense-fluid  approach  is  that  the  residence  time  is  very  short  because  of 
rapid  mixing  with  the  hot  reaction  products,  which  quickly  pushes  the  mixture  well  above  the  vapor  dome. 
The  present  work  considers  the  numerical  and  physical  modeling  issues  that  arise  when  this  “single-fluid” 
approach  is  used  for  representing  transcritical  flows  under  both  non-reacting  and  reacting  conditions. 

Other  more  comphensive  modeling  options  are  certainly  possible  but  they  are  much  more  involved.  A 
homogenous  two  fluid  model  considers  the  vapor  state  and  liquid  state  to  be  independent  species.  The 
species  are  assumed  to  be  uniformly  mixed  and  at  the  same  temperature.  There  is  no  interface  tracking  or 
surface  tension  in  this  model.2  In  the  case  when  surface  tension  becomes  important,  Eulerian-Lagrangian 
approaches  that  prescribe  spray  distributions  and  breakup  can  be  used.3  When  empirically  based  atomization 
and  breakup  models  are  insufficient,  a  complete  interface  tracking  method  like  the  volume  of  fluid  approach 
can  be  used.  In  such  an  approach  the  liquid  gas  interface  is  tracked.4  We  note  that  an  overwhelming  majority 
of  the  computational  studies  have  similarly  focused  on  purely  supercritical  flows,  which  serve  to  validate  the 
use  of  different  real  gas  equations  of  state,  but  fail  to  address  the  more  complex  transcritical  regime. 

Over  the  past  decade,  a  substantial  amount  of  work  has  been  directed  at  studying  the  behavior  of  jets  in 
transcritical  environments.  Chehroudi  et  al.  collected  visual  images  and  jet  spreading  rates  for  jet  injection 
into  like  and  unlike  ambient  fluids  for  subcritical  and  supercritical  conditions.5,6  Leyva  et  al.  have  looked 
at  supercritical  nitrogen  injection  into  a  nitrogen  environment  under  both  unforced  and  acoustically  forced 
conditions.7  Roy  et  al.  performed  experimental  work  on  the  injection  of  supercritical  jets  into  subcritical 
environments.  Depending  on  the  thermodynamic  state  they  were  able  to  produce  results  where  surface 
tension  became  important  and  droplets  were  formed.  Other  conditions  showed  no  droplet  formation  and 
behavior  similar  to  supercritical/supercritical  injection.8  This  indicates  that  this  is  an  extremely  complex 
operating  state  without  clearly  defined  flow  regimes. 

Previous  modeling  work  by  Okong’o  et  al.  focused  on  a  supercritical  binary  mixing  in  a  mixing  layer. 
Using  direct  numerical  simulation  (DNS)  they  found  that  the  strong  density  gradient  hindered  entrainment. 
Results  also  showed  that  the  turbulent  diffusion  of  species  was  more  important  than  momentum.9  Further 
work  by  Masi  et  al.  on  turbulent  mixing  layers  developed  effective  locally  varying  turbulent  Schmidt  and 
Prandtl  numbers  which  were  observed  to  have  negative  values.10  These  two  parameters  can  be  used  control 
the  turbulent  diffusion  of  mass  and  energy  and  have  typically  been  assumed  to  be  constant  positive  values 
in  low-pressure  environments. 

In  a  comprehensive  review  article,  Bellan  cites  two  supercritical  jet  simulations  and  notes  that  a  real 
gas  equation  of  state  along  with  unsteady  simulations  are  the  absolute  minimum  requirements  for  capturing 
supercritical  behavior.11  More  recently,  Zong  et  al.  simulated  supercritical  N2  injected  into  gaseous  N2  and 
found  good  agreement  in  the  jet  spreading  angle  with  experimental  data  for  two  density  ratios.  We  note 
that  an  overwhelming  majority  of  the  computational  studies  have  similarly  focused  on  purely  supercritical 
flows,  which  serve  to  validate  the  use  of  different  real  gas  equations  of  state,  but  fail  to  address  the  more 
complex  transcritical  regime.  Oefelein  simulated  combustion  between  transcritical  oxygen  and  supercritical 
hydrogen.  It  was  found  that  intense  property  gradients  were  present  that  approached  contact  discontinuities 
which  can  be  difficult  to  model.12  Recently,  Dahms  and  Oefelein  have  proposed  a  theoretical  framework  for 
assessing  when  a  propellant  mixture  behaves  as  a  supercritical  dense  fluid  and  when  transcritical  phenomena 
needs  to  be  considered.13  However,  modeling  of  the  transcritical  regime  remains  a  challenge. 

As  noted  above,  we  evaluate  the  use  of  a  single-fluid  model  based  on  the  Peng  Robinson  equation  of 
state  under  supercritical  and  transcritical  conditions.  Although  our  ultimate  interest  is  in  combustion,  we 
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consider  both  non-reacting  and  reacting  flows  in  order  to  understand  the  underlying  modeling  challenges. 
Under  subcritical  conditions,  the  intrinsic  discontinuity  present  in  the  equation  of  state  constitutes  a  physical 
and  numerical  inconsistency  with  the  continuum  assumption.  For  this  reason,  our  objective  is  to  modify 
the  equation  of  state  to  remove  the  singularity  and  relax  it  to  a  continuous  variation  in  thermodynamic 
phase  space.  Accordingly,  the  paper  is  organized  as  follows.  In  the  next  section,  we  summarize  the  Peng 
Robinson  equation  of  state  and  present  several  modifications  and  smoothing  treatments  designed  to  keep 
the  system  of  equations  well-behaved.  In  addition,  we  consider  different  mixing  rules  in  order  to  better 
control  the  mixture  properties.  We  show  that  each  of  these  methods  have  strengths  and  weaknesses.  We 
consider  results  for  both  non-reacting  and  reacting  jets.  The  non-reacting  test  cases  are  drawn  from  the 
experimental  configurations  of  Cheroudi  et  al.  that  were  mentioned  erlier,  while  the  reacting  configuration 
corresponds  to  operating  conditions  and  geometries  relevant  to  liquid  rocket  engines.  Based  on  these  results, 
we  provide  several  conclusions  about  the  relative  efficacy  of  employing  single-fluid  models  for  supercritical 
and  transcritical  problems. 


II.  Equation  of  State  Model  Adjustments 


A.  Peng  Robinson  Equation  of  State 

A  large  number  of  cubic  equations  of  state  are  available.  A  summary  of  the  numerous  different  models 
is  available  in  Poling  et  al.14  The  cubic  equation  of  state  used  in  the  present  work  is  the  Peng- Robinson 
equation  of  state.1  This  is  a  three-parameter  cubic  equation  of  state  which  takes  the  form, 


RT  dm 

Vm-bm  +  2brnVrn  ~ 


(1) 


The  values  of  am  and  bm  can  be  found  through  a  series  of  mixing  rules.  First,  defining  the  values  for  an 
individual  species  i, 


Ki  =  0.37464  +  1.54226^  -  0.26992^ 
R2T^ 

a*  =  0.457235 - ^ 

Pc,i 

bi  =  0.077796 

Pc,i 


2 


(2) 

(3) 

(4) 


The  parameters  am  and  bm  are  calculated  through  the  following  mixing  rules  using  the  values  obtained  for 
the  individual  species  in  the  equations  above, 


N  N 

& vn  ^  ^  Xj Xj  yjdjdj  j  (5) 

j= 1 i=l 

N 

brn  =  ^2,Xibi  (6) 

2  =  1 


Here  k'-  is  an  interaction  parameter  which  can  be  used  to  adjust  the  mixing  rule  for  pairs  of  species,  in  most 
cases  it  is  unity.  Equation  (1)  can  alternatively  be  written  as  a  cubic  equation,  the  roots  of  which  are  the 
compressibility  Z, 


Z3  -  (1  -  B)Z2  +  (A-2B-  W2)Z  -  ( AB  -  B2  —  B3)  =  0 


(7) 


where  the  compressibility  is, 


pW 

Jet 


(8) 


In  the  case  of  multiple  positive  real  roots,  the  largest  value  of  Z  is  typically  selected.  The  parameters  A  and 
B  depend  on  the  mixture  and  are  related  to  dm  and  bm  through  the  following  definitions, 


A-a  P 

171  R2T2 

(9) 

B  =  bm -P 

RT 

(10) 
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In  addition  to  the  above  equations,  departure  functions  are  used  to  calculate  all  other  thermodynamic 
properties.14  Figure  1  shows  the  density  of  nitrogen  calculated  using  the  Peng  Robinson  equation  of  state 
for  several  reduced  pressures.  The  model  captures  the  high  density  at  very  low  temperatures  and  the 
supercritical  behavior  at  high  reduced  pressures.  The  model  does  not  however  provide  a  solution  to  the 
thermodynamic  state  inside  the  dome.  Reduced  pressures  of  0.43,  0.63  and  0.83  show  sharp  transitions 
from  the  high  density  to  the  low-density  regimes.  The  implication  of  the  single  fluid  model  is  that  near  the 
transition  point  when  the  temperature  shifts  slightly  the  density  (and  other  properties)  will  immediately 
jump  in  value.  This  presents  a  numerical  challenge  because  small  changes  in  temperature  (e.g.,  less  than 
0.1  K)  can  result  in  order  of  magnitude  changes  in  density  and  specific  heat.  Above  a  reduced  pressure  of  1.03 
the  transition  between  the  high-density  region  and  low-density  region  is  smooth,  but  still  exhibits  non-ideal 
gas  behavior.  A  similar  plot  can  be  generated  for  a  multi-component  mixture.  The  challenge  with  a  mixture 
is  that  the  transition  region  is  also  a  function  of  composition  in  addition  to  the  pressure  and  temperature. 
In  a  reacting  flow  computation  avoiding  this  problematic  transition  region  is  impossible  due  to  the  wide 
varying  properties  of  the  mixture. 


100  150  200  250  300 

Temperature,  K 

Figure  1:  Density  over  the  temperature  range  of  interest.  Each  curve  corresponds  to  a  different  reduced 
pressure. 


B.  Compressibility  Smoothing 

In  addition  to  the  sharp  jump  of  density  in  the  vapor  dome  at  subcritical  conditions  (reduced  pressure  lower 
than  unity),  the  use  of  a  cubic  equation  of  state  in  the  context  of  a  single  fluid  model  leads  to  a  discontinuity  in 
other  thermodynamic  properties.  Enthalpy  and  sound  speed  are  no  longer  continuous  and  constant  pressure 
specific  heat  tends  to  infinity.  This  results  in  numerical  stiffness  that  cannot  be  handled  properly  by  most 
of  numerical  methods  currently  used  in  computational  fluid  dynamics  in  the  space  propulsion  community. 

The  initial  proposed  solution  to  improve  the  numerical  stability  while  remaining  with  a  single  fluid  model 
is  to  artificially  smooth  the  transition  region  through  the  dome  region.  The  argument  being  that  in  practice 
the  exact  behavior  inside  this  region  for  injector  type  simulations  is  restricted  to  a  very  small  region  in  the 
flowfield.  Creating  a  method  to  traverse  this  region  ensures  that  the  simulation  progresses  in  a  stable  fashion. 
In  the  dome  region  there  will  be  three  real  roots  to  the  cubic  equation  (Eq.  7).  Let  Zi  be  the  ordered  roots 
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to  Eq.  7  ordered  such  that, 


Zx>  Z2>  Z3  (11) 

Figure  2  shows  a  plot  of  the  three  roots  in  the  region  of  interest  for  nitrogen  at  a  reduced  pressure  of  0.63. 
The  density  is  also  shown  for  comparison.  The  compressibility  used  to  define  the  density  is  always  Z\ ,  which 
is  ma x{Z^}.  In  the  region  near  the  discontinuity  there  are  three  roots  present,  the  largest  root  corresponds 
to  the  vapor  phase  density  and  the  smallest  root  corresponds  to  the  liquid  phase  density.  The  middle  root 
is  without  physical  meaning.  The  temperature  range  over  which  the  three  roots  are  present  is  bounded  by 
Tl  on  the  low  side  and  Th  on  the  high  side. 


Figure  2:  Roots  to  equation  7  along  with  density  for  nitrogen  at  a  reduced  pressure  of  0.63. 


In  the  current  formulation,  the  compressibility  is  artificially  smoothed  through  the  transition  region  using 
the  following  definition  of  a  “smooth”  compressibility  Zs, 

Zs  =  Zi-  —  ~  Z:i  [1  +  tanh  G]  (12) 

where  the  parameter  G  is, 

G  =  ag  (P7*Zl  --Z,  +  (Z2  -  Z3)\  (13) 

f  IT  Pmax  J 

Figure  3  shows  the  smooth  root  and  the  corresponding  density.  The  thermodynamic  derivatives  of  Z  are 
also  needed,  and  it  was  found  that  not  applying  the  appropriate  smoothing  to  derivatives  of  Z  resulted  in 
undesirable  behavior  of  other  thermodynamic  quantities,  especially  at  the  interface  where  Zs  switches  to  Z\. 
The  derivative  of  the  smooth  root  is, 


G'  = 


Z[~-  (Z[  -  Z'3)  (1  +  tanh  G)  +  (Z1  -  Z3)  G' 


cosh2  G 


Zf2-Z[ 


1 


Pm  ax 


Pmax  T  1  1  T  P  n 


(14) 

(15) 


The  smooth  density  works  very  well  for  a  single  species.  In  the  case  of  a  mixture  it  was  found  that  away 
from  the  dome  region  the  smoothing  could  alter  the  solution.  Figure  4  shows  a  mixture  of  70%  C12H26, 
10%  C02,  10%  02  and  10%  H2  for  a  variety  of  pressures.  Notice  how  in  the  two  highest  pressure  cases 
at  around  1300  K  there  is  a  spike  in  the  density.  A  closer  examination  of  solution  to  the  cubic  equation 
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Figure  3:  Roots  to  equation  7  along  with  the  smooth  root  and  the  smoothed  density  for  nitrogen  at  a  reduced 
pressure  of  0.63. 


in  this  region  reveals  the  presence  of  three  real  roots  with  at  least  one  negative  root.  In  order  to  improve 
the  robustness  of  the  smoothing  and  to  allow  it  to  be  applied  universally  the  criterion  was  modified  to  only 
apply  the  smoothing  where  all  three  real  roots  were  found  to  be  positive.  Making  this  change  restores  the 
correct  behavior  of  the  high-pressure  solutions  at  1300  K,  as  shown  in  the  figure. 

The  smoothing  formulation  introduces  two  new  parameters,  pm ax  and  ag.  pm ax  controls  the  location  of 
the  smoothing  and  ag  controls  the  transition  curvature.  The  smoothing  is  fairly  insensitive  to  the  parameter 
ag  and  a  value  of  30  seems  to  provide  adequate  curvature  at  the  transition  points.  Figure  5  shows  the  effect 
of  pm ax  on  pure  nitrogen  at  a  reduced  pressure  of  0.23.  By  tuning  pmax  the  location  of  the  transition  is 
shifted.  The  smoothing  alters  other  properties  including  the  specific  heat  and  sound  speed.  The  expected 
behavior  for  the  specific  heat  is  to  have  a  very  sharp  maximum  at  the  critical  temperature;  this  is  instead 
shifted  with  the  transition  region,  and  the  maximum  specific  heat  shifts  from  being  a  near  singular  value 
to  a  broad  peak  with  an  overall  maximum  less  than  the  unaltered  solution.  The  sound  speed  should  reach 
a  local  minimum  at  the  critical  temperature.  In  the  case  of  the  smoothing  the  minimum  sound  speed  is 
further  reduced  and  the  transition  occurs  over  a  larger  temperature  range.  The  sound  speed  also  exhibits 
some  discontinuity  where  the  compressibility  shifts  from  Zs  to  Z\  at  the  high  temperature  point. 

The  difficulty  with  this  smoothing  is  that  it  introduces  unwanted  changes  to  other  thermodynamic  prop¬ 
erties  like  the  specific  heat  and  sound  speed.  Additionally,  it  is  difficult  to  choose  a  single  value  of  pmax  that 
is  suitable  for  a  wide  range  of  pressures.  This  is  further  complicated  in  the  multi-species  cases  where  the 
transition  region  is  a  function  of  temperature,  pressure,  and  species  composition.  To  address  this  several 
additional  smoothing  approaches  were  considered.  The  following  sections  briefly  describe  each  approach. 
Following  the  descriptions  a  comparison  of  each  of  the  methods  is  given.  The  problematic  transition  region 
where  three  real  positive  roots  exist  can  be  bounded  by  a  lower  temperature  Tl  and  an  upper  temperature 
Th.  This  range  AT  is  the  entire  range  there  the  three  root  system  exists  and  was  shown  in  Figure  2.  The 
solution  to  a  general  cubic  equation, 

x 3  +  ax 2  +  bx  +  c  =  0  (16) 

the  discriminant  determines  the  types  of  roots,  real  or  complex,  and  is  defined  as, 

A  =  i  (2«3  -  9ab  +  27)2  -  -h  +  3b)3  (17) 
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70%  Cl2 10%  CO;  10%  O;  10%  H;0 


70%  ClzHze  10%  CO;  10%  02  10%  H;0 


Temperature.  K 


Temperature.  K 


(a)  Smoothing  applied  any  time  three  real  roots  were  (b)  Smoothing  applied  only  when  three  positive  real 
found.  roots  were  found. 


Figure  4:  A  negative  real  root  can  cause  the  smoothing  operation  to  fail  outside  the  dome  region  and  alter 
the  original  solution. 


The  presence  of  three  real  roots  occurs  when  the  discriminant  A  is  less  than  zero.  To  identify  the  temperature 
range  under  which  the  three  real  roots  system  is  active  the  sign  of  the  discriminant  is  examined.  While 
varying  temperature  for  a  fixed  pressure  and  mass  composition  the  first  temperature  where  the  sign  of  the 
discriminant  changes  corresponds  to  Tl  and  where  the  discriminant  changes  sign  again  the  temperature 
is  Th .  Figure  6  shows  the  identification  of  the  low  and  high  temperatures  based  on  the  discriminant  for 
nitrogen  at  a  reduced  pressure  of  0.63. 

1.  Modified  Hyperbolic  Tangent  Function 

The  first  alternative  method  still  uses  a  hyperbolic  tangent  to  form  the  transition.  At  Tl  and  Th  the 
compressibility  is  a  smooth  function.  In  an  effort  to  preserve  the  smooth  property  at  these  points  which 
were  shown  to  be  problematic  in  the  prior  definition  the  smooth  root  is  defined  in  terms  of  only  the  maximum 
Z\  and  the  minimum  Z3  compressibility.  The  smooth  compressibility  is  defined  as, 


^  [Z\  (1  -  tanh 6S)  +  Z3  (1  +  tanh#s)] 

(18) 

where 

T  —  Tc 

(19) 

0S  =  4 

Th  -  T 

and  the  center  temperature  Tc  is, 

m  _Tl+Th 
c  ~  2 

(20) 

2.  Mirrored  Z2 


The  roots  already  form  a  continuous  path  between  the  minimum  and  maximum  compressibility  but  it  is  not 
a  one-to-one  function.  To  exploit  the  continuous  path,  the  root  Z2  must  be  mirrored  so  that  Z  is  a  function. 
The  smooth  root  is  defined  as  a  piece- wise  function  of  the  three  roots, 


Z. a  —  < 


Zi 


3A  T 
4 

AT 


<  T  <  Th 


mirrored  Z2  ^  <  T  <  ^ 


Us 


TL  <T  < 


AT 


(21) 
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(a)  Density.  (b)  Sound  Speed. 


(c)  Specific  Heat. 

Figure  5:  The  value  of  pmax  is  changed  for  pure  N2  at  a  reduced  pressure  of  0.23.  In  each  case  the  black 
dashed  line  represents  the  original  Peng  Robinson  solution. 


8  of  20 


American  Institute  of  Aeronautics  and  Astronautics 


Figure  6:  Variation  of  the  discriminant  of  the  cubic  equation. 


3.  Cubic  Polynomial  of  Temperature 

Another  straightforward  way  of  approximating  the  compressibility  is  by  obtaining  a  cubic  fit  between  the 
compressibility  at  the  lower  and  higher  temperatures.  To  fit  the  coefficients  of  the  cubic  function  Z(Tl ) 
and  Z(Th )  are  used.  Variations  in  temperature  produce  large  changes  in  Z  compared  with  pressure,  and  to 
account  for  this  and  are  used  as  the  additional  equations  to  compute  the  coefficients  of  the 

cubic  equation.  The  smooth  root  is  then, 

Zs  =  O-i  +  02  T  +  (23  T2  +  (24  T3  (22) 

with  the  parameters  obtained  using  the  boundary  conditions  described  above. 

4-  Using  a  Cosine  Function 

The  analytic  solution  of  the  cubic  equation  in  the  region  of  interest  requires  determining  a  quantity  0.  This 
0  can  be  shown  to  vary  within  the  interval  [0,  |] .  The  earlier  use  of  hyperbolic  tangent  requires  the  domain 
to  be  [-4,4]  for  obtaining  [-1,1]  range.  Instead,  the  theta  variation  in  this  case  can  be  exploited  to  formulate 
a  smooth  transition  as, 


Zs  =  t[Z3(  1  +  cos  39)  +  Zi  (1  —  cos  36>)]  (23) 

The  results  of  all  of  these  approaches  along  with  the  original  formulation  have  strengths  and  weaknesses 
in  terms  of  specific  properties.  The  difficulty  is  that  there  is  not  a  clear  “best  choice”  in  smoothing  the 
roots.  Figure  7  shows  a  comparison  of  all  proposed  methods  in  determining  the  density,  sound  speed,  and 
specific  heat  for  C12H26  at  a  pressure  of  1.8  MPa  for  a  temperature  range  that  spans  the  region  of  interest. 
Qualitatively  the  density  is  similar  in  each  case,  taking  a  slightly  different  path  for  each.  The  modified 
tanh  and  cosine  approach  have  sharper  transitions,  especially  at  The  sound  speed  variations  show 
very  different  behavior  from  each  other  and  from  the  original  Peng  Robinson  solution.  The  cosine,  cubic 
polynomial,  and  modified  tanh  approach  cause  the  sound  speed  to  increase  in  the  transition  region  unlike 
the  physical  behavior.  The  original  tanh  and  mirrored  Z2  method  have  a  lower  sound  speed  that  is  shifted 
away  from  the  critical  temperature.  While  the  mirrored  Z2  method  showed  promise  in  the  density  and  sound 
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speed  it  generates  infinite  values  for  the  specific  heat.  That  is  because  where  Z\  transitions  to  Z2  and  where 
Z2  transitions  to  Z3  the  temperature  derivative  of  compressibility  is  infinite.  This  effectively  rules  out  this 
method.  All  of  the  other  smoothing  methods  have  a  large  band  of  high  specific  heats  shifted  away  from 
the  Peng  Robinson  value.  Aside  from  eliminating  the  mirrored  Z 2  method  there  is  no  clear  best  smoothing 
approach. 


Figure  7:  The  value  of  pmax  is  changed  for  pure  N2  at  a  reduced  pressure  of  0.23.  In  each  case  the  black 
dashed  line  represents  the  original  Peng  Robinson  solution. 


C.  Combining  Peng  Robinson  with  Amagat’s  Mixing  Rule 

The  complex  mixing  rules  associated  with  the  principle  of  corresponding  states  results  in  problematic  con¬ 
ditions  when  small  amounts  of  water  or  carbon  dioxide  are  added  to  cold  fuel.  This  is  a  situation  that  will 
regularly  occur  in  liquid  rocket  injectors  where  cold  fuel  slowly  begins  to  burn  and  small  amounts  of  water 
or  carbon  dioxide  are  present  at  lower  temperatures.  To  help  mitigate  this  a  modified  mixing  rule  based 
on  Amagat’s  law  has  been  implemented.  In  this  approach,  each  species  is  treated  individually  as  a  real 
gas.  This  means  that  instead  of  a  global  compressibility,  each  species  will  have  its  own  compressibility.  In 
addition,  departure  functions  for  other  properties  are  computed  on  a  per-species  basis  as  opposed  to  one 
departure  function  for  the  mixture.  The  following  mixing  rules  are  consistent  with  Amagat’s  law  of  partial 
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volumes.  The  mixture  density  is, 


/  N  \  1 

II 

SM 

cii 

(24) 

where  p,  is  the  partial  density,  defined  as, 

pWi 

Pl  RTZi 

(25) 

Other  properties  like  enthalpy  use  the  following  mixing  rule, 

N 

h  =  J2(hi+  hP)  Yi 

(26) 

2=1 


Where  hi  is  the  ideal  gas  enthalpy  of  species  i  and  h f*  is  the  enthalpy  departure  function  for  species  i.  Other 
properties  like  specific  heat,  entropy,  and  Gibbs  energy  are  evaluated  in  the  same  manner.  The  mixtures 
compressibility  can  be  defined  as, 

N 

Z  =  YJZiXi  (27) 

2=1 

This  approach  yields  a  mixture  compressibility  that  is  similar  to  that  found  using  the  standard  Peng  Robin¬ 
son  model.15  This  model  offers  another  potential  advantage  in  controlling  previously  observed  problematic 
behavior.  Each  species  can  be  smoothed  individually  potentially  removing  the  challenge  in  defining  a  uni¬ 
versal  smoothing  parameter  applicable  to  a  disparate  mixture.  Figure  8  shows  a  comparison  of  a  mixture  of 
Ci2H26  with  either  water  or  carbon  dioxide  at  300  K  for  3  MPa  and  10  MPa.  For  water,  there  is  very  little 
difference  in  the  computed  densities  for  the  two  mixing  rules.  Carbon  dioxide  shows  a  significant  departure 
at  the  lower  pressure.  This  is  because  the  range  of  C12H26  mass  fractions  spans  the  dome  region. 


(a)  C12H26  and  C02  at  300  K.  (b)  C12H26  and  H20  at  300  K. 

Figure  8:  The  value  of  pmax  is  changed  for  pure  N2  at  a  reduced  pressure  of  0.23.  In  each  case  the  black 
dashed  line  represents  the  original  Peng  Robinson  solution. 
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D.  Constant  Compressibility  for  Reacting  Shear  Layers 

Similar  to  the  prior  model  in  this  case  we  have  assumed  that  the  real  gas  effects  are  limited  to  a  very  narrow 
region  of  the  flowfield.  The  key  issue  with  using  an  ideal  gas  model  is  that  the  density  of  the  fuel  (or  oxidizer) 
will  be  too  low.  Since  the  mass  flow  rate  is  typically  specified  as  the  inlet  boundary  condition  using  the 
wrong  inlet  density  will  result  in  velocities  that  are  incorrect.  In  this  model  a  fixed  value  of  Z  is  specified 
according  to, 

Z=^~  (28) 

Pig 

Here  the  real  gas  density,  prg ,  is  computed  using  the  inlet  boundary  conditions  from  the  Peng  Robinon 
equation  of  state  or  a  database  like  REFPROP.  The  mixing  rule  used  in  equation  24  is  used.  This  model  is 
more  limited  in  that  no  real  gas  effects  in  the  properties  other  than  density  are  accounted  for,  and  the  value 
of  Z  is  only  accurate  at  the  inlet  conditions,  as  the  fuel  or  oxidizer  deviate  from  this  operating  point  the 
approximation  will  worsen.  In  the  configurations  of  interest  however  the  fuel  occupies  only  a  small  region 
and  is  typically  quickly  consumed.  Figure  9  shows  a  comparison  of  the  compressibility  computed  with  Peng 
Robinson  and  the  constant  compressibility  supplied  at  300  K  and  600  K. 


Figure  9:  Comparision  of  the  compressibility  computed  with  Peng  Robinson  and  the  constant  compressibility 
supplied  at  300  K  and  600  K. 

Based  on  the  figure  it  is  clear  that  this  approximation  is  only  valid  when  the  fuel  remains  close  to  the 
injection  temperature  where  the  compressibility  was  defined.  Because  the  ideal  gas  curve  does  not  match 
the  real  gas  curve  a  simple  multiplication  of  the  compressibility  has  limited  validity  outside  the  region  where 
it  was  defined.  In  a  liquid  rocket  injector,  we  expect  that  the  fuel  is  quickly  burned  and  is  present  only  in  a 
small  portion  of  the  flowfield  thus  an  such  approximation  may  be  acceptable.  Importantly,  the  combustion 
products  are  a  a  higher  temperature  wherein  the  ideal  gas  assumption  is  usually  valid. 
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III.  Non-reacting  Results 


The  simulations  for  the  non-reacting  case  follow  the  work  of  Chehroudi  et  al.  who  performed  a  series  of 
experiments  primarily  using  nitrogen.  Cold  N2  (90  —  110  K)  is  injected  into  a  chamber  filled  with  warm  N2 
(300  K).  This  was  done  for  a  series  of  pressures  where  the  reduced  pressure  ranged  from  0.23  to  2.74.  This 
spans  the  subcritical  to  supercritical  pressure  states.  In  doing  so,  they  recorded  vastly  different  injection 
behavior.  At  the  low  reduced  pressure  the  jet  is  liquid  like  with  surface  instabilities.  Above  a  reduced  pressure 
of  unity  there  is  no  breakup  into  droplets,  instead  there  is  a  mixing  layer  with  large  density  uniformities.  At 
the  highest  reduced  pressure  the  jet  resembles  a  single  phase  turbulent  jet.6 

A  representative  two-dimensional  planar  geometry  was  created  to  examine  the  applicability  of  the  smooth¬ 
ing  to  capture  the  subcritical  to  supercritical  pressure  states.  A  5.5  mm  by  0.4064  mm  domain  uses  a  uniform 
256  x  512  Cartesian  grid.  The  code  is  fifth  order  accurate  in  space  and  a  third  order  explicit  RungeKutta 
scheme  is  used  for  the  time  discretization.  The  simulations  are  run  without  a  turbulence  model.  Nitrogen  at 
90  K  is  injected  through  a  height  of  0.127  mm  at  a  velocity  of  12.5  m/s.  The  remainder  of  the  left  boundary 
flows  nitrogen  at  300  K  with  a  velocity  of  0.5  m/s.  Simulations  were  attempted  using  the  tanh  smoothing 
with  the  parameter  ag  =  30  and  pmax  =  4.0  for  reduced  pressures  of  0.43,  0.63,  0.83,  1.03,  1.26,  1.64,  2.03, 
2.44,  and  2.74.  All  simulations  with  a  reduced  pressure  of  1.03  or  less  failed  to  generate  a  stable  solution 
that  lasted  through  a  single  flow  through  time.  Simulations  with  a  reduced  pressure  of  1.26  and  above  were 
stable. 

Figure  10  shows  density  contours  for  the  reduced  pressure  of  1.64  and  Figure  11  shows  density  contours  for 
the  reduced  pressure  of  2.74.  Both  simulations  show  a  clear  demarcation  between  the  warm  and  cold  nitrogen. 
The  interface  between  the  two  fluid  streams  is  characterized  by  a  Kelvin  Helmholtz  type  instability.  In  the 
case  of  the  reduced  pressure  of  1.64  the  roll-up  is  larger  and  starts  closer  to  the  inlet  boundary  compared 
with  the  reduced  pressure  of  2.74.  Experimentally  for  the  jet  it  was  reported  that  above  a  reduced  pressure 
of  2.04  there  was  a  reduction  in  the  core  length  and  by  2.74  it  visually  looked  like  gas-gas  mixing.  The  shear 
layer  test  case  makes  it  difficult  to  draw  conclusions  about  the  core  length  but  the  results  do  visually  look 
like  gas-gas  mixing. 
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Figure  10:  Density  plots  for  a  reduced  pressure  of  1.64. 
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Figure  11:  Density  plots  for  a  reduced  pressure  of  2.74. 
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IV.  Reacting  Results 


A  two-dimensional  shear  layer  representative  of  a  C12H26  injector  is  used  to  test  the  constant  compress¬ 
ibility  ideal  gas  model.  For  this  simulation,  a  code  which  is  second  order  accurate  in  time  and  space  is  used. 
The  mesh  contains  approximately  250,000  cells  and  has  uniform  spacing.  The  fuel  is  C12H26  and  is  injected 
at  600  K  with  a  mass  flow  rate  of  74.085  kg/s  •  m,  the  oxidizer  is  pure  02  and  is  injected  at  611  K  with  a 
mass  flow  rate  of  145.3567  kg/s  •  m.  The  mean  pressure  is  3.846  MPa.  For  the  fuel,  these  conditions  are 
supercritical  in  pressure  but  subcritical  in  temperature.  For  the  oxidizer,  the  conditions  are  subcritical  in 
pressure  and  supercritical  in  temperature.  A  temperature  of  600  K  was  chosen  for  the  fuel  because  under 
simplified  conditions  without  acoustic  waves  this  case  can  be  successfully  run  using  the  Peng  Robinson  model 
with  tanh  smoothing  (ag  =  30  and  pmax  =  4.0).  The  actual  configuration  is  subjected  to  hydrodynamic 
and  combustion  instabilities  and  thus  generates  wide  ranging  pressure  fluctuations  that  are  not  currently 
numerically  stable  with  the  Peng  Robinson  model.  A  simplified  two  step  reaction  mechanism  is  used. 

Results  from  the  two  simulations  are  showing  in  Figures  12-14.  The  figures  show  the  density,  fuel  mass 
fraction,  and  temperature.  From  the  density  plot  it  is  clear  that  the  fixed  compressibility  has  difficulty 
maintaining  a  near  uniform  incoming  density.  Looking  at  different  time  snapshots  of  the  density  (not 
shown)  shows  that  the  density  varies  unlike  the  real  gas  approach  which  shows  the  same  density  at  both 
times.  This  is  likely  because  of  pressure  oscillations  and  the  fact  that  the  ideal  gas  equation  of  state  with  a 
fixed  compressibility  does  not  follow  the  same  density  curve  as  the  real  gas  as  the  temperature  and  pressure 
change  from  the  evaluation  point  of  Z.  Both  simulations  start  burning  at  slightly  different  times  so  it  is 
difficult  to  compare  direct  snapshots. 


x,  m 

(a)  Peng  Robinson. 


Density,  kg/m3 


20  60  100  140  ISO  220  260  300  340  360  420  460 


x,  m 

(b)  Constant  Z. 

Figure  12:  Density  plots  for  the  reacting  shearlayer  test  case. 

The  fuel  mass  fraction  along  with  the  density  identify  differences  in  how  the  shear  layer  interface  behaves 
for  the  two  models.  One  difference  is  the  sharpness  of  the  interface  between  the  fuel  and  oxygen.  The  real 
gas  shows  a  steep  gradient  between  the  two  while  the  density  gradients  are  smoothed  out  in  the  other  case. 
This  is  a  result  of  the  equation  of  state  and  should  be  further  considered  in  evaluating  the  applicability  of 
this  model  and  how  it  affects  the  mixing  of  the  propellants.  There  is  also  large  scale  roll  up  present  in  the 
fixed  compressibility  case  that  is  not  observed  in  the  real  gas  simulation.  Another  difference  is  that  as  the 
fuel  heats  up  (downstream)  the  density  is  significantly  lower  for  the  real  gas  equation  of  state.  This  is  to 
be  expected  based  on  Figure  9,  which  showed  that  the  density  would  be  over  predicted  based  on  evaluating 
the  compressibility  at  600  K.  A  piecewise  construction  of  Z  at  multiple  temperatures  may  be  able  to  better 
capture  the  behavior 

Despite  these  differences  the  temperature  plot  shows  a  qualitatively  similar  combustion  with  a  thin  flame 
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along  the  interface  between  the  fuel  and  oxygen.  The  flame  is  attached  in  both  cases  and  is  anchored  to  the 
splitting  plate.  This  model  is  an  attractive  engineering  model,  especially  when  the  computed  value  of  Z  is 
close  to  unity.  But  there  are  indications  that  under  some  conditions  it  may  introduce  significant  differences 
from  the  real  gas  equation  of  state.  A  full  understanding  of  how  the  fuel  behaves  and  the  length  of  time 
it  remains  in  the  combustor  are  key  questions  to  answer  before  deciding  if  this  model  is  applicable  to  a 
particular  configuration. 


x,  m 


(a)  Peng  Robinson. 


x,  m 

(b)  Constant  Z. 


Figure  13:  C12H26  mass  fraction  plots  for  the  reacting  shearlayer  test  case. 
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Figure  14:  Temperature  plots  for  the  reacting  shearlayer  test  case. 
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V.  Conclusions 


A  series  of  smoothing  options  were  presented  to  enable  the  Peng  Robinson  equation  of  state  to  perform 
without  introducing  numerically  destablizing  jumps  in  properties  in  the  transcritical  regime.  The  motivation 
was  to  enable  a  single  fluid  model  for  applications  in  liquid  rocket  injectors  where  the  transcritical  regime  is 
restricted  to  a  small  local  region  inside  the  injector  cup  where  the  fuel  and  oxidizer  initially  begin  to  burn. 
The  smoothing  proved  to  be  difficult  and  focused  on  the  computation  of  the  compressibility  in  the  vapor 
dome.  While  it  was  simple  to  generate  a  density  that  had  unique  values  for  temperature  and  pressure,  other 
thermodynamic  quantities  were  problematic.  Both  the  constant  pressure  specific  heat  and  the  sound  speed 
have  unique  behavior  near  the  critical  point.  The  sound  speed  reaches  a  local  minimum  and  the  specific  heat 
takes  on  extremely  large  values  over  a  very  narrow  range.  The  smoothing  was  unable  to  fully  capture  that 
behavior  despite  multiple  different  approaches  to  generate  a  smooth  compressibility. 

A  non-reacting  nitrogen  shear  layer  was  tested  for  a  range  of  reduced  pressures  using  hyperbolic  tangent 
smoothing  that  provided  the  best  overall  behavior.  It  was  found  that  only  simulations  above  the  vapor  dome 
(i.e.,  in  the  compressed  liquid  regime)  were  able  to  run  successfully.  Simulations  which  were  sub-critical  in 
pressure  were  numerically  unstable.  This  may  indicate  that  the  single  fluid  approach  to  modeling  this  regime 
is  inadequate  and  a  Eulerian-Lagrangian  or  VOF  approach  is  needed.  It  may  also  be  a  result  of  inadequate 
numerical  stability  to  handle  the  extremely  large  thermodynamic  gradients  that  are  present.  This  will  be 
considered  in  future  work.  Another  issue  is  that  the  smoothing  introduces  inconsistencies,  the  value  of  Z 
is  not  consistent  with  am  and  bm  in  the  smoothing  region.  Since  all  three  of  these  are  used  to  compute 
other  thermodynamic  quantities  it  may  be  necessary  to  consistently  recompute  am  and  bm.  This  is  not 
straightforward  and  will  also  be  considered  in  future  studies. 

While  the  non-reacting  case  is  interesting,  a  major  motivation  of  this  work  is  the  configuration  with 
combustion.  In  some  regards,  this  appears  to  be  a  simpler  problem.  As  the  flowfield  reacts,  the  increase  in 
temperature  naturally  moves  it  away  from  the  vapor  dome.  The  challenge  for  the  reacting  flow  case  is  the 
introduction  of  water  and  carbon  dioxide,  which  because  of  their  high  critical  pressure,  tend  to  shift  mixtures 
toward  the  vapor  dome.  This  can  create  local  regions  with  vastly  different  properties,  including  sound  speeds 
so  low  that  the  flow  becomes  locally  supersonic.  Two  additional  equation  of  state  modifications  were  identified 
for  this  type  of  simulation.  A  modified  mixing  rule  which  used  Amagat’s  law  and  a  fixed  compressibility. 
Both  of  these  approaches  exploit  the  fact  that  the  real  gas  behavior  is  restricted  to  the  incoming  propellants 
and,  for  a  short  portion  of  time  over  a  limited  area,  when  the  propellants  mix.  The  majority  of  the  flowfield 
where  combustion  has  taken  place  will  behave  more  like  an  ideal  gas.  The  fixed  compressibility  case  was 
compared  with  a  smooth  Peng  Robinson  case  in  a  reacting  shear  layer  configuration.  The  results  showed 
qualitatively  similar  burning.  There  were  differences  in  the  computed  density  the  Peng  Robinson  equation  of 
state  produced  a  near  uniform  density  of  the  incoming  fuel  while  the  fixed  compressibility  showed  variability 
because  of  local  pressure  fluctuations.  Further  work  is  needed  to  understand  how  to  better  control  this, 
possibly  by  moving  to  the  modified  mixing  rule. 
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